rm(list = ls())

# Parameters
C.cost <- 0.5
Overline.v <- 5
Shape.1 <- 2 
Shape.2 <- 2
Mean <- 2
Sd <- 2

# Functions
W.0 <- function(v){
  out <- 1-pnorm(q = v, mean = 0, sd = 5)
  return(out) 
}

W.1 <- function(t, w1.max = 1, w1.slope = 0.05){
  out <- w1.max - t*w1.slope
  return(out)
}

find.equilibrium <- function(t, 
                             cost = C.cost,
                             mean = Mean, sd = Sd,
                             overline.v = Overline.v,
                             shape1 = Shape.1, shape2 = Shape.2){
  L <- length(t)
  lower <- rep(NA, L)
  equ.threshold <- rep(NA, L)
  f <- function(v, t){
    q <- (1-pnorm(overline.v, mean = mean, sd = sd)) / (1-pnorm(v, mean = mean, sd = sd))
    admission <- pbeta(q = q, shape1 = shape1, shape2 = shape2)
    adj.stakes <- cost / (W.1(t = t) - W.0(v))
    out <- admission - adj.stakes
    return(out)
  }
  for(i in 1:L){
    lower[i] <- 5*qnorm(1-W.1(t[i])+cost)
    upper <- overline.v
    equ.threshold[i] <- uniroot(f, 
                                lower = lower[i], 
                                upper = upper, 
                                t = t[i])$root
  }
  data <- data.frame(Equ.threshold = equ.threshold, 
                     Lower = lower)
  return(data)
}

utility <- function(t, 
                    coeff = 0.02, exponent = 1, 
                    mean = Mean, sd = Sd,
                    shape1 = Shape.1, shape2 = Shape.2, 
                    cost = C.cost, overline.v = Overline.v){
  L <- length(t)
  application <- rep(NA, L)
  posterior <- rep(NA, L)
  admission <- rep(NA, L)
  v.hat <- rep(NA, L)
  utility <- rep(NA, L)
  for(i in 1:L){
    v.hat[i] <- find.equilibrium(t = t[i],  
                                 mean = mean, sd = sd, 
                                 shape1 = shape1, shape2 = shape2,
                                 cost = cost, overline.v = overline.v)$Equ.threshold
    application[i] <- 1 - pnorm(q = v.hat[i], mean = mean, sd = sd)
    posterior[i] <- (1 - pnorm(overline.v, mean = mean, sd = sd)) / (1 - pnorm(v.hat[i], mean = mean, sd = sd))
    admission[i] <- pbeta(posterior[i], shape1 = shape1, shape2 = shape2)
    utility[i] <- 1 - application[i]*admission[i] - (coeff*t[i]^exponent)
  }
  data <- data.frame(Utility = utility, 
                     Threshold = v.hat, 
                     Application = application, 
                     Admission = admission, 
                     Posterior = posterior)
  return(data)
}

# Figure

t.vary <- seq(0, 0.9, length.out = 1000)

pdf('./figures/fig_e3.pdf', width = 14, height = 7)
par(mfrow = c(1, 2))
plot(t.vary, utility(t = t.vary, coeff = 0.001,
                     exponent = 17,
                     shape2 = 4.85)$Utility, 
     type = 'l', lwd = 5,
     yaxt = 'n', 
     xlab = 'Asylum Law', ylab = '',
     cex.lab = 1.35, col = 'darkgray', ylim = c(0.8728471, 0.8729298))
title(ylab = 'Utility Politician', cex.lab = 1.35, line = 1.3)
lines(t.vary, utility(t = t.vary, coeff = 0.001,
                      exponent = 12,
                      shape2 = 4.85)$Utility, 
      col = 'gray', lwd = 5)
lines(t.vary, utility(t = t.vary, coeff = 0,
                      shape2 = 4.85)$Utility, 
      col = 'black', lwd = 5, lty = 2)
legend(x = 0, y = 0.87292, 
       lwd = c(3, 3, 3), lty = c(2, 1, 1), 
       bty = 'n',
       col = c('black', 'gray', 'darkgray'), 
       legend = c('No costs', 
                  expression(0.001%.%t^12), 
                  expression(0.001%.%t^17)), 
       cex = 1.5)

plot(t.vary, utility(t = t.vary, coeff = 0.001,
                     exponent = 17,
                     shape2 = 6)$Utility, 
     type = 'l', lwd = 5, 
     yaxt = 'n', col = 'darkgray',
     xlab = 'Asylum Law', ylab = '',
     cex.lab = 1.35, ylim = c(0.8505323, 0.8511307))
title(ylab = 'Utility Politician', cex.lab = 1.35, line = 1.3)
lines(t.vary, utility(t = t.vary, coeff = 0.001,
                      exponent = 12,
                      shape2 = 6)$Utility, 
      col = 'gray', lwd = 5)
lines(t.vary, utility(t = t.vary, coeff = 0,
                      shape2 = 6)$Utility, 
      col = 'black', lwd = 5, lty = 2)
legend(x = 0, y = 0.85105, 
       lwd = c(3, 3, 3), lty = c(2, 1, 1),
       bty = 'n',
       col = c('black', 'gray', 'darkgray'), 
       legend = c('No costs', 
                  expression(0.001%.%t^12), 
                  expression(0.001%.%t^17)),
       cex = 1.5)
dev.off()
